home *** CD-ROM | disk | FTP | other *** search
/ Practical Algorithms for Image Analysis / Practical Algorithms for Image Analysis.iso / CH_5.7 / EH_SEG / EH_SEG.C < prev    next >
C/C++ Source or Header  |  1999-09-11  |  14KB  |  578 lines
  1. /*
  2.  * eh_seg.c
  3.  *
  4.  * Practical Algorithms for Image Analysis
  5.  *
  6.  * Copyright (c) 1997, 1998, 1999 MLMSoftwareGroup, LLC
  7.  */
  8.  
  9. /*
  10.  * ** E(xtract) H(istograms)_SEG
  11.  * **
  12.  * ** routine to extract desired set of parameters from a data file of 
  13.  * ** type .SEG whose name is entered as a command line argument and 
  14.  * ** construct a histogram; write output file of type *.hdt to be read
  15.  * ** by routine hf.c
  16.  * **
  17.  */
  18. #include "eh_seg.h"
  19.  
  20. #define MAX_REC_SIZE    1200
  21.  
  22. #define BIN_NUMBER    10           /* no bins in histogram */
  23. #define    LEN_HIST    "\nhistogram for segment length"
  24. #define    ANGLE_HIST    "\nhistogram for segment turn angles"
  25.  
  26. #define    ON         1
  27. #define    OFF         0
  28.  
  29. #undef    SHOW_ANG
  30. #undef    SHOW_SORT
  31. #define    DBG_ANG
  32.  
  33.  
  34. /*
  35.  * ** global variables
  36.  */
  37. extern char *optarg;
  38. extern int optind, opterr;
  39.  
  40.  
  41. static char *HIST_HEADER;
  42. int N_BINS = BIN_NUMBER;
  43. int MIN_SET = 0;
  44. int MAX_SET = 0;
  45. int WRITE_FILE = 0;
  46. int SHOW_INPUT = 0;
  47. int A_HIST = 0;
  48. int L_HIST = 0;
  49.  
  50. /*
  51.  * usage of routine
  52.  */
  53. void
  54. usage (char *progname)
  55. {
  56.   progname = last_bs (progname);
  57.   printf ("USAGE: %s infile [-s] [-a] [-l] [-n n] [-i f] [-f f] [-w] [-L]\n", progname);
  58.   printf ("\n%s extracts histogram data from infile (SEG), generated by\n", progname);
  59.   printf ("fitpolyg, containing line segment data\n\n");
  60.   printf ("ARGUMENTS:\n");
  61.   printf ("  infile: input filename (SEG)\n\n");
  62.   printf ("OPTIONS:\n");
  63.   printf ("    -s: show input data\n\n");
  64.   printf ("    construct histogram for:\n");
  65.   printf ("    -a: turn angles between consecutive segments\n");
  66.   printf ("    -l: histogram of segment lengths\n");
  67.   printf ("  -n n: set number of bins to n (default: 10)\n");
  68.   printf ("  -i f: set initial value to (float)f (default: MIN)\n");
  69.   printf ("  -f f: set final value to (float)f (default: MAX)\n");
  70.   printf ("    -w: write output file of type .hdt (hist data)\n");
  71.   exit (1);
  72. }
  73.  
  74.  
  75.  
  76. /*
  77.  * ** comparison function for qsort():
  78.  * ** sort array of tuples in order of increassing x-values
  79.  */
  80. int 
  81. compare (t1, t2)
  82.      float *t1, *t2;
  83. {
  84.   return ((int) SIGN (*t1 - *t2));
  85. }
  86.  
  87.  
  88. /*
  89.  * ** determine mean of input data set
  90.  */
  91. double 
  92. find_mean (float *data, int n)
  93. {
  94.   int i;
  95.   double mean = 0.0;
  96.  
  97.   for (i = 0; i < n; i++)
  98.     mean += *(data + i);
  99.  
  100.   return (mean / (double) n);
  101. }
  102.  
  103. /*
  104.  * ** determine standard deviation of input data set
  105.  */
  106. double 
  107. find_sigma (float *data, int n, double mean)
  108. {
  109.   int i;
  110.   double xi, sigma = 0.0;
  111.  
  112.   for (i = 0; i < n; i++) {
  113.     xi = (double) (*(data + i));
  114.     sigma += (xi - mean) * (xi - mean);
  115.   }
  116.   sigma /= (double) (n - 1);
  117.  
  118.   return (sqrt (sigma));
  119. }
  120.  
  121.  
  122. /*
  123.  * ** construct histogram of input data
  124.  */
  125. void 
  126. construct_hist (int n, float *data, float *hist, double bw, double data_base)
  127. {
  128.   int i, ibin;
  129.  
  130.   for (i = 0; i < n; i++) {
  131.     if (*(data + i) != 0) {
  132.       ibin = 0;
  133.       while (*(data + i) >= data_base + ibin * bw)
  134.         ibin++;
  135.       *(hist + ibin - 1) += 1;
  136.     }
  137.   }
  138. }
  139.  
  140.  
  141.  
  142. /*
  143.  * ** read first line in data file to determine size of record
  144.  */
  145. int 
  146. seg_size (FILE * file)
  147. {
  148.   int retval;
  149.   int n_rec;
  150.  
  151.   if (((retval = fscanf (file, "%d", &n_rec)) == 0) || (retval == EOF)) {
  152.     printf ("wrong input file format!\n");
  153.     exit (1);
  154.   }
  155.   return (n_rec);
  156. }
  157.  
  158.  
  159. #define    LOLIM          0.01
  160. #define    UPLIM        200.0
  161.  
  162. /*
  163.  * ** decode given slope, tan PHI, and return signed angle PHI (in rad);
  164.  * ** sign convention: CCW <--> pos, counting from pos. x-direction
  165.  */
  166. double 
  167. slp_to_angle (double m, double delx, double dely)
  168. {
  169.   double phi = -5.0;
  170.  
  171. #ifdef DEBUG
  172.   printf ("SLP_TO_ANGLE:...m = %lf, delx = %lf, dely = %lf\n",
  173.           m, delx, dely);
  174. #endif
  175.  
  176.   if (fabs (m) < LOLIM) {       /* horizontal segm */
  177.     if (SIGN (delx) > 0)
  178.       phi = 0.0;
  179.     else if (SIGN (delx) < 0)
  180.       phi = PI;
  181.     else
  182.       phi = -7.0;
  183.   }
  184.   else if (fabs (m) > UPLIM) {  /* vertical segm */
  185.     if (SIGN (dely) < 0)
  186.       phi = 0.5 * PI;
  187.     else if (SIGN (dely) > 0)
  188.       phi = 1.5 * PI;
  189.     else
  190.       phi = -8.0;
  191.   }
  192.   else {
  193.     if (SIGN (delx) > 0) {
  194.       if (SIGN (dely) > 0)
  195.         phi = atan (m);
  196.       if (SIGN (dely) < 0)
  197.         phi = 2.0 * PI + atan (m);
  198.     }
  199.     else if (SIGN (delx) < 0) {
  200.       if (SIGN (dely) > 0)
  201.         phi = PI + atan (m);
  202.       if (SIGN (dely) < 0)
  203.         phi = PI + atan (m);
  204.     }
  205.     else
  206.       phi = -9.0;
  207.   }
  208.  
  209.   if (fabs (phi) > 2.0 * PI) {
  210.     printf ("\nSLP_TO_AN:...ill defined condition encountered:");
  211.     printf (" m = %lf, phi = %lf\n", m, phi);
  212.     phi = 0.0;
  213.     return (phi);
  214.   }
  215.  
  216. #ifdef DEBUG
  217.   printf ("                value to be returned: %8.4lf\n", phi);
  218. #endif
  219.  
  220.   return (phi);
  221. }
  222.  
  223.  
  224. /*
  225.  * ** read segment data from formatted data file (generally of type .seg)
  226.  * ** written by pcctoseg.c; first column contains segment index: must not
  227.  * ** exceed 9999!!
  228.  */
  229. void 
  230. get_seg_data (fp, n_segm, segm)
  231.      FILE *fp;
  232.      int n_segm;
  233.      struct Segm *segm;
  234. {
  235.   int i;
  236.   int retval;
  237.  
  238.   for (i = 0; i < n_segm; i++) {
  239.     retval = fscanf (fp, "%4d %3d %3d %3d %3d %f %3d",
  240.                      &((segm + i)->segm_ind),
  241.                      &((segm + i)->ptO.x), &((segm + i)->ptO.y),
  242.                      &((segm + i)->ptF.x), &((segm + i)->ptF.y),
  243.                      &((segm + i)->slope), &((segm + i)->line_ind));
  244.   }
  245. }
  246.  
  247.  
  248.  
  249.  
  250. void 
  251. main (argc, argv)
  252.      int argc;
  253.      char *argv[];
  254. {
  255.   int i_arg;
  256.   int i, ich, is;
  257.   int n_segm;
  258.   static char wbuf[20];
  259.   static char *inp_suffix =
  260.   {".seg"};                     /* default inp file suffix */
  261.   static char *seg_type =
  262.   {".seg"};                     /* suffix for .sgl inp file */
  263.   static char *wsuffix =
  264.   {".hdt"};                     /* suffix for output file name */
  265.   char *rbuf;
  266.  
  267.   int n_parms = 3;
  268.   double ndel_x, ndel_y, del_x, del_y;
  269.   struct Segm *inSegm;          /* sample input */
  270.  
  271.   int cur_ln_ind, next_ln_ind;
  272.   float cur_slp, next_slp;
  273.   double cur_th, next_th, del_th;
  274.  
  275.  
  276.   float bin_width, min = (float) -1.0, max = (float) -1.0;
  277.   float *data, *hist;
  278.   double mean, std_dev;
  279.   float one = (float) 1.0;
  280.  
  281.   FILE *fpIn, *fpOut;
  282.  
  283. /* 
  284.  * ** cmd line options:
  285.  */
  286.   static char *optstring = "slan:i:f:wL";
  287.  
  288.  
  289. /*
  290.  * ** parse command line
  291.  */
  292.   optind = 2;
  293.   opterr = ON;                  /* give error messages */
  294.  
  295.  
  296.   if (argc < 2)
  297.     usage (argv[0]);
  298.  
  299.   while ((i_arg = getopt (argc, argv, optstring)) != EOF) {
  300.     switch (i_arg) {
  301.     case 's':
  302.       printf ("\n...option %c: ", i_arg);
  303.       printf ("display input data\n");
  304.       SHOW_INPUT = ON;
  305.       break;
  306.     case 'a':
  307.       printf ("\n...option %c: ", i_arg);
  308.       printf ("hist of turn angles betw consec segm\n");
  309.       A_HIST = ON;
  310.       break;
  311.     case 'l':
  312.       printf ("\n...option %c: ", i_arg);
  313.       printf ("hist of segm lengths\n");
  314.       L_HIST = ON;
  315.       break;
  316.     case 'n':
  317.       printf ("\n...option %c: ", i_arg);
  318.       printf ("set number of bins to %d\n",
  319.               N_BINS = atoi (optarg));
  320.       break;
  321.     case 'i':
  322.       printf ("\n...option %c: ", i_arg);
  323.       printf ("set min (init value) to %f\n",
  324.               min = (float) atof (optarg));
  325.       MIN_SET = 1;
  326.       break;
  327.     case 'f':
  328.       printf ("\n...option %c: ", i_arg);
  329.       printf ("set max (final value) to %f\n",
  330.               max = (float) atof (optarg));
  331.       MAX_SET = 1;
  332.       break;
  333.     case 'w':
  334.       printf ("\n...option %c: ", i_arg);
  335.       printf ("write output to file %s\n", wbuf);
  336.       WRITE_FILE = 1;
  337.       break;
  338.     case 'L':
  339.       print_sos_lic ();
  340.       exit (0);
  341.     default:
  342.       printf ("\n...unknown cmd line argument\n");
  343.       exit (1);
  344.       break;
  345.     }
  346.   }
  347.   printf ("read sgl data from file %s\n", rbuf = argv[1]);
  348.   /* open input file */
  349.   if ((fpIn = fopen (rbuf, "r")) == NULL) {
  350.     printf ("\n...cannot not open input file\n");
  351.     exit (1);
  352.   }
  353.  
  354.   /* initialize size variables */
  355.   if ((n_segm = seg_size (fpIn)) > MAX_REC_SIZE) {
  356.     printf ("\n...record size of %d exceeds limit of %d\n",
  357.             n_segm, MAX_REC_SIZE);
  358.     exit (1);
  359.   }
  360.  
  361.   /* construct output file name */
  362.   ich = 0;
  363.   while ((*(wbuf + ich) = *(rbuf + ich)) != '.')
  364.     ich++;
  365.   for (is = 0; is < 4; is++)
  366.     *(wbuf + (ich) + is) = *(wsuffix + is);
  367.   /* strip input file suffix */
  368.   for (is = 0; is < 4; is++)
  369.     *(inp_suffix + is) = *(rbuf + (ich) + is);
  370.   if (strcmp (inp_suffix, seg_type) != 0) {
  371.     printf ("\n...unknown input file type encountered:");
  372.     printf (" default: %s\n", seg_type);
  373.     fclose (fpIn);
  374.     exit (1);
  375.   }
  376.   if (WRITE_FILE) {
  377.     if ((fpOut = fopen (wbuf, "w")) == NULL) {
  378.       printf ("\n...cannot open %s for writing\n", wbuf);
  379.       exit (1);
  380.     }
  381.   }
  382.   if ((A_HIST == 0) && (L_HIST == 0)) {
  383.     printf ("\n...must specify option -a or -l!!\n");
  384.     exit (1);
  385.   }
  386.   if (A_HIST == 1) {
  387.     if (min > max) {
  388.       printf ("\n...delected min exceeds max!!\n");
  389.       exit (1);
  390.     }
  391.   }
  392.   printf ("...number of segments in file %s: %d\n", rbuf, n_segm);
  393.  
  394.  
  395. /*
  396.  * ** allocate memory
  397.  */
  398.   if ((hist = (float *) calloc ((size_t) N_BINS, sizeof (float))) == NULL) {
  399.     printf ("\nEH_SEG: mem alloc for hist failed\n");
  400.     exit (1);
  401.   }
  402.   if ((inSegm = (struct Segm *) calloc ((size_t) n_segm, sizeof (struct Segm))) == NULL) {
  403.     printf ("\nEH_SEG: mem alloc for inSegm failed\n");
  404.     exit (1);
  405.   }
  406.   if ((data = (float *) calloc ((size_t) n_segm, sizeof (float))) == NULL) {
  407.     printf ("\nEH_SEG: mem alloc for data failed\n");
  408.     exit (1);
  409.   }
  410.  
  411. /*
  412.  * ** fetch data
  413.  */
  414.   get_seg_data (fpIn, n_segm, inSegm);
  415.  
  416. /*
  417.  * ** write output
  418.  */
  419.   if (SHOW_INPUT == ON) {
  420.     printf ("\n...input data:\n");
  421.     for (i = 0; i < n_segm; i++) {
  422.       printf ("%4d %3d %3d %3d %3d %f8.2 %4d\n",
  423.               (inSegm + i)->segm_ind,
  424.               (inSegm + i)->ptO.x, (inSegm + i)->ptO.y,
  425.               (inSegm + i)->ptF.x, (inSegm + i)->ptF.y,
  426.               (inSegm + i)->slope, (inSegm + i)->line_ind);
  427.     }
  428.   }
  429.  
  430. /*
  431.  * ** evaluate statistics and construct histogram
  432.  */
  433.   if (L_HIST == ON) {
  434.     printf ("\n...evaluating histogram of segment lengths...\n");
  435.     for (i = 0; i < n_segm; i++) {
  436.       del_x = (double) ((inSegm + i)->ptF.x - (inSegm + i)->ptO.x);
  437.       del_y = (double) ((inSegm + i)->ptF.y - (inSegm + i)->ptO.y);
  438.       *(data + i) = (float) sqrt (del_x * del_x + del_y * del_y);
  439.     }
  440.   }
  441.   else if (A_HIST == ON) {
  442.     printf ("\n...evaluating histogram of turn angles...\n");
  443.     i = 0;
  444.     do {
  445.  
  446.       cur_ln_ind = (inSegm + i)->line_ind;
  447. #ifdef SHOW_ANG
  448.       printf ("\nline_index: %4d\n", cur_ln_ind);
  449. #endif
  450.  
  451.       while ((next_ln_ind = (inSegm + i + 1)->line_ind) == cur_ln_ind) {
  452. #ifdef SHOW_ANG
  453.         printf ("...segment pair:%4d %4d ",
  454.                 (inSegm + i)->segm_ind, (inSegm + i + 1)->segm_ind);
  455. #endif
  456.         cur_slp = (inSegm + i)->slope;
  457.         del_x = (double) ((inSegm + i)->ptF.x - (inSegm + i)->ptO.x);
  458.         del_y = (double) ((inSegm + i)->ptF.y - (inSegm + i)->ptO.y);
  459.         next_slp = (inSegm + i + 1)->slope;
  460.         ndel_x = (double) ((inSegm + i + 1)->ptF.x - (inSegm + i + 1)->ptO.x);
  461.         ndel_y = (double) ((inSegm + i + 1)->ptF.y - (inSegm + i + 1)->ptO.y);
  462.  
  463.         /* angle conv: CCW <--> positive */
  464.         if (cur_slp * next_slp == -1.0) {
  465.           printf ("\n...segments perpendicular...\n");
  466.           if (cur_slp > 0.0) {
  467.             if (del_x > 0.0) {
  468.               if (ndel_x > 0.0)
  469.                 *(data + i) = (float) 90.0;
  470.               if (ndel_x < 0.0)
  471.                 *(data + i) = (float) -90.0;
  472.             }
  473.             else if (del_x < 0.0) {
  474.               if (ndel_x > 0.0)
  475.                 *(data + i) = (float) -90.0;
  476.               if (ndel_x < 0.0)
  477.                 *(data + i) = (float) 90.0;
  478.             }
  479.           }
  480.           else if (cur_slp < 0.0) {
  481.             if (del_x > 0.0) {
  482.               if (ndel_x > 0.0)
  483.                 *(data + i) = (float) -90.0;
  484.               if (ndel_x < 0.0)
  485.                 *(data + i) = (float) 90.0;
  486.             }
  487.             else if (del_x < 0.0) {
  488.               if (ndel_x > 0.0)
  489.                 *(data + i) = (float) 90.0;
  490.               if (ndel_x < 0.0)
  491.                 *(data + i) = (float) -90.0;
  492.             }
  493.           }
  494.         }
  495.         else {
  496.           cur_th = slp_to_angle (cur_slp, del_x, del_y);
  497.           next_th = slp_to_angle (next_slp, ndel_x, ndel_y);
  498.           del_th = next_th - cur_th;
  499. #ifdef SHOW_ANG
  500.           printf (" th: %5.2lf nth: %5.2lf d_th: %5.2lf",
  501.                   cur_th, next_th, del_th);
  502. #endif
  503.  
  504.           if (del_th > PI)
  505.             del_th = -(2.0 * PI - del_th);
  506.           if (del_th < -PI)
  507.             del_th = 2.0 * PI + del_th;
  508.  
  509.           *(data + i) = (float) ((del_th * 180.0) / PI);
  510.  
  511. #ifdef SHOW_ANG
  512.           printf (" d_th(deg): %5.2f\n", *(data + i));
  513. #endif
  514.         }
  515.         cur_ln_ind = next_ln_ind;
  516.         i++;
  517.       }
  518.       i++;
  519.     } while (i < n_segm);
  520.   }
  521.  
  522.  
  523.   printf ("\n...sorting and constructing histogram\n");
  524. #if defined(LINUX)
  525.   qsort (data, n_segm, sizeof (float), (__compar_fn_t) compare);
  526. #else
  527.   qsort (data, n_segm, sizeof (float), compare);
  528. #endif
  529.  
  530. #ifdef SHOW_SORT
  531.   printf ("\n...sorted data:\n");
  532.   for (i_sgl = 0; i_sgl < n_segm; i_sgl++)
  533.     printf (" data[%d] = %f\n", i_sgl, *(data + i_sgl));
  534. #endif
  535.  
  536.   if (MIN_SET == 0)
  537.     min = *(data + 0);
  538.   if (MAX_SET == 0)
  539.     max = *(data + n_segm - 1);
  540.  
  541.   mean = find_mean (data, n_segm);
  542.   std_dev = find_sigma (data, n_segm, mean);
  543.  
  544.   bin_width = (max - min) / ((float) (N_BINS - 1));
  545.   construct_hist (n_segm, data, hist, bin_width, min);
  546.  
  547.   printf ("\nhistogram:\n");
  548.   printf ("......mean: %lf\n", mean);
  549.   printf ("...std_dev: %lf\n", std_dev);
  550.   printf ("...min: %f, bin_width: %f, max: %f\n", min, bin_width, max);
  551.   for (i = 0; i < N_BINS; i++) {
  552.     printf (" %6.2f", *(hist + i));
  553.     if ((i + 1) % 10 == 0)
  554.       printf ("\n");
  555.   }
  556.   fclose (fpIn);
  557.  
  558. /*
  559.  * ** write output file of type .hdt (see also: .adt, .pdt)
  560.  */
  561.   if (WRITE_FILE == 1) {
  562.     fprintf (fpOut, "%d\n", N_BINS);
  563.  
  564.     for (i = 0; i < n_parms; i++)
  565.       fprintf (fpOut, "%f\n", one);
  566.  
  567.     for (i = 0; i < N_BINS; i++)
  568.       fprintf (fpOut, "%f %f\n", (float) i, *(hist + i));
  569.  
  570.     fclose (fpOut);
  571.   }
  572.  
  573.  
  574.   free (hist);
  575.   free (inSegm);
  576.   free (data);
  577. }
  578.